Understanding the emergence of the boson peak in molecular glasses

A common feature of glasses is the “boson peak”, observed as an excess in the heat capacity over the crystal or as an additional peak in the terahertz vibrational spectrum. The microscopic origins of this peak are not well understood; the emergence of locally ordered structures has been put forward as a possible candidate. Here, we show that depolarised Raman scattering in liquids consisting of highly symmetric molecules can be used to isolate the boson peak, allowing its detailed observation from the liquid into the glass. The boson peak in the vibrational spectrum matches the excess heat capacity. As the boson peak intensifies on cooling, wide-angle x-ray scattering shows the simultaneous appearance of a pre-peak due to molecular clusters consisting of circa 20 molecules. Atomistic molecular dynamics simulations indicate that these are caused by over-coordinated molecules. These findings represent an essential step toward our understanding of the physics of vitrification.


Summary
In this paper is examined the origin of the boson peak using a prototypical molecular glass former, tetrabutyl orthosilicate (TBOS). The experimental methods consist of optical Kerr effect (OKE) spectroscopy as well as wide-and small-angle x-ray scattering (SAXS,WAXS) measurements, complemented by other several other measurements (such as calorimetry, NMR, FTIR, Raman, viscometry) and simulations (MD and DFT). The main results relate to an excess temperature dependent intensity in the <3 THz region of the Raman spectrum, probed by OKE. The spectrum is analysed in terms of inter-and intramolecular components; the former consists of a high-frequency band at ~2THz attributed to acoustic phonons and a second low-frequency band ~0.9 THz attributed to the boson peak. The SAXS/WAXS measurements indicate the existence of a pre-peak at ~0.234 /Å which is attributed to the formation of molecular clusters corresponding to distances of ~26.9 Å. The MD simulations reproduce well both the OKE and SAXS/WAXS results, whereas Voronoi analysis indicates that below Tg = 124 K (estimated by calorimetry) emerge clusters with BCC local order. In conclusion, the observed band attributed to the boson peak is associated with the formation of over-coordinated clusters which can prevent crystallization. Assessment I find the paper well written, with clarity of presentation and rich as well as convincing results. I liked the described scheme for using probing the boson peak by the OKE approach by selectivity through symmetry of the selected system (here near tetrahedral), which allows to exclude rotational contributions (confirmed by DFT). I also find convincing that this assignment is confirmed by low-temperature calorimetry which shows an excess heat capacity at 6K. The results are of very high quality, including a nice comparison between the OKE measurements with WAXS measurements. The additional measurements are almost overwhelming, although useful for strengthening the interpretation of transient cluster formation while at the same time showing that the sample remains monomeric, as confirmed by NMR and DFT. In the results is also included MD simulations that allow an interesting insight into the molecular coordination of the clusters, which are mainly icosahedral (at any given temperature) and truncated octahedral (appearing at T_g) . Understanding the origin of the boson peak can be a of great interest for many in the field and gives a promising insight for understanding the link of collective oscillations due to "cage rattling" modes which manifest in the supercooled region, which can inhibit crystallization. I therefore would support the manuscript's publication in Nature Communications after the following comments have been considered. Comments 1. I was wondering about the interpretation of the black line shown in Fig. 1b, which shows a pronounced temperature dependence as indicated in Suppl. Fig. 6. I missed assignment of this peak in the text, even though this is the range with the most pronounced T-dependence (at least in the log-x representation).
2. I have a few comments related to the pre-peak indicated in the WAXS at ~0.234 /Å: 2a. As far as I understand the WAXS and SAXS spectra are "stitched together" just below ~0.2 /Å, so it would be good to specify the exact q-range in the methods section.
2b. I am a bit concerned about the influence of the background either from the capillary or the direct beam. It is mentioned that plastic capillaries were used, so it would be useful to plot the SAXS/WAXS background from the capillaries (without any sample) to ensure that there is no matching peak at this momentum transfer range.
2c. Suppl. Fig. 8 needs some additional care. It would be helpful to make the figure labels visible, include proper legends and also re-bin (using larger q-bins) the SAXS data to improve the readability of the data in the q-range above 0.01/Å.
3. In connection to the previous comment, I am also wondering about the pre-peak observed S(q) simulations results shown in Fig. 4b and Suppl. Fig. 14: 3a. While indeed as it is pointed out in the paper, the pre-peak is reproduced in the Si-Si S(q), similar peaks in the region q<1/Å with even stronger temperature dependence are shown in the S(q)'s of C-O, C-Si and also in the O-O. Do these temperature-dependent signals reflect changes intermolecular or intramolecular distances? I.e. it should be specified whether in the calculation is included atoms within the same molecule or only intermolecular distances.
3b. It might be helpful to plot all these results in Suppl. Fig. 14 with logarithmic x-axis (as in Fig. 4b).
3c. The peak at q~0.2 /Å is attributed to clusters with size ~2π/q = 26.9 Å. Would then one not expect a similar peak in the Si-Si PCF at this range (Suppl. Fig.15)? Or is this range limited by the simulation box size? 4. Concerning the Voronoi analysis: 4a. It would be helpful to include some additional explanations concerning the formalism used in the corresponding indices and how one can derive the corresponding symmetry.
4b. It would be interesting to include a schematic of the Si-SI network and polyhedral derived from the simulations to facilitate the interpretation.
4c. The analysis indicates that network of the Si atoms exhibits icosahedral coordination at any given temperature, whereas the observed truncated octahedra emerge near the T_g. One of the main conclusions of the paper is however that the observed transient clusters can inhibit crystallization. According to this interpretation, wouldn't one expect that the cluster formation should emerge in the supercooled regime?

Reviewer #2 (Remarks to the Author):
The work by González-Jiménez studied the dynamical properties, i.e. boson peak, of a molecular glass by combining experiments and simulations. A link to the local structure is suggested. However, this paper is more technical relevant with lots of assumptions involved in the analysis, so it is not convincing to me. Therefore, I cannot recommend its publication in Nat. Commun., a more specialized journal may be more suitable. See some detailed comments below. 1) The introduction should be carefully re-organized to target the problem to be solved in this paper. For example, the main results are related to boson peak and glass structure, but the introduction mainly focused on different relaxation modes, like beta relaxation. Not many details about the boson peak are shown.
2) In the introduction, "Below the glass-transition, secondary or β-relaxation processes occur due to faster dynamic processes" may not be necessarily true. I believe different relaxation modes should coexist always in a thermodynamic state, but their observation will depend on the experimental time scale.
3) There is a contradiction about the origin of the secondary relaxations in the introduction, whether intramolecular or intermolecular. 4) The detailed data analysis from the experimental spectra depends on the assumption that they will be simplified for some molecules with ordinary symmetry, tetrahedral and octahedral. This has not been rationalized. Also, in the glass state, there are enormous types of these local structures with different distortions, how would this fact influence the results? 5) The spectrum fitting with four functions of many free parameters seems highly possible to bring over-fitting. This is very critical for the overall interpretation of the results. 6) The structure analysis by pair-correlation function and the Voronoi tessellation does not account for the orientation information in the local structures, which should be important in this molecular glass. 7) The direct comparison of the simulations by a simplified model with the experiments is not convincing. The consistency of the partial structure factor only shows the overall averaged information similarly, not necessarily mean the model it's correct. 8) At elevated temperatures, anharmonicity is important, which should have been into consideration in the spectrum analysis.

Reviewer #3 (Remarks to the Author):
The manuscript of González-Jiménez et al. entitled "Understanding the emergence of the boson peak in molecular glasses" aims at unveiling the boson peak in a highly symmetric glass former using femtosecond optical Kerr-effect spectroscopy as the central characterization technique, at different temperatures. They then use other experimental techniques (WAXS, SAXS, Raman) and simulation to characterize this boson peak. It is because the characterized molecule, TBOS, is symmetric that the authors can tackle this longstanding problem of the origin of the boson peak. Two intermolecular bands in the OKE spectrum have been mainly analyzed. In Figure 1 d, the amplitudes of both bands are shown versus the temperature, and in Figure 2 a, their ratio is displayed. Through theses analyses, the authors assign the low frequency band as the boson peak, and the high frequency band as the acoustic phonon band. Thanks to this assignment, they surmised that the boson peak is associated with the formation of clusters of TBOS molecules undergoing collective oscillation. WAXS experiments lead to a distance of 26.9 Angstrom for clusters, and from MD simulation, 25.1 Angstroms. Conclusions drawn from the OKE experiments can have been made because of the symmetric nature of the TBOS molecule. Regarding the over-coordinated molecules in the clusters, deduced from MD, the authors wrote: "the MD simulations are strongly suggestive that these structures are clusters of over-coordinated TBOS molecules". For an article in Nature Communications, the authors must be more convincing. For instance, as the molecular behavior can be probed, the over-coordination can be proved. Other issues. p.3: It is written "The calculated infrared spectrum in the Si-O-C stretch region (around 1,100 cm-1) matches the experimental one (see Supplementary Figure 3)". Only the 1000-1200 cm-1 spectrum should be considered. A great part of the simulated spectrum does not fit the experimental one. The authors should indicate why the fit works only for the region of interest. p.4: the Supplementary Figure 9 does not correspond to the liquid density as outlined in the text. Indicate what are q3 and q5. Below 110 K: the value seems to be 110 K. It is written that "the TBOS molecules reduce the number of gauche defects in the alkoxide side chains". Such a behavior should be observed in MD. p.5: The authors consider the two coordination shells in the RDF, computing 2 and 12 molecules respectively. According to these spectra, these numbers should change with the temperature. Some relevant results can be obtained. p.9: MD simulation part. When the authors consider random positions within a cubic box: what is the random algorithm used? It is necessary to specify that the periodic boundary conditions have been applied. Simulations are performed in the NPT ensemble, meaning that a thermostat and a barostat are needed simultaneously. So, why the authors wrote: "to sample the canonical ensemble"? The authors wrote: "In this regime, the self-diffusion is sufficient to sample a reasonable portion of the configurational space for the liquid phase". One point in the configurational space corresponds to 87,552 coordinates. What do the authors mean with this reasonable portion?
Finally, how can the study presented in this paper be correlated with the recent paper in Nature Physics on the origin of the boson peak? " Reply to reviewer comments Reviewer #1 (Remarks to the Author):

Summary
In this paper is examined the origin of the boson peak using a prototypical molecular glass former, tetrabutyl orthosilicate (TBOS). The experimental methods consist of optical Kerr effect (OKE) spectroscopy as well as wide-and small-angle x-ray scattering (SAXS,WAXS) measurements, complemented by other several other measurements (such as calorimetry, NMR, FTIR, Raman, viscometry) and simulations (MD and DFT). The main results relate to an excess temperature dependent intensity in the <N THz region of the Raman spectrum, probed by OKE. The spectrum is analysed in terms of inter-and intramolecular components; the former consists of a high-frequency band at ~TTHz attributed to acoustic phonons and a second low-frequency band ~U.V THz attributed to the boson peak. The SAXS/WAXS measurements indicate the existence of a pre-peak at ~U.TNX /Å which is attributed to the formation of molecular clusters corresponding to distances of ~TZ.V Å. The MD simulations reproduce well both the OKE and SAXS/WAXS results, whereas Voronoi analysis indicates that below Tg = ]TX K (estimated by calorimetry) emerge clusters with BCC local order. In conclusion, the observed band attributed to the boson peak is associated with the formation of over-coordinated clusters which can prevent crystallization. Assessment I find the paper well written, with clarity of presentation and rich as well as convincing results. I liked the described scheme for using probing the boson peak by the OKE approach by selectivity through symmetry of the selected system (here near tetrahedral), which allows to exclude rotational contributions (confirmed by DFT). I also find convincing that this assignment is confirmed by low-temperature calorimetry which shows an excess heat capacity at ZK. The results are of very high quality, including a nice comparison between the OKE measurements with WAXS measurements. The additional measurements are almost overwhelming, although useful for strengthening the interpretation of transient cluster formation while at the same time showing that the sample remains monomeric, as confirmed by NMR and DFT. In the results is also included MD simulations that allow an interesting insight into the molecular coordination of the clusters, which are mainly icosahedral (at any given temperature) and truncated octahedral (appearing at T_g) . Understanding the origin of the boson peak can be a of great interest for many in the field and gives a promising insight for understanding the link of collective oscillations due to "cage rattling" modes which manifest in the supercooled region, which can inhibit crystallization.
We thank the reviewer for the very supportive comments.
I therefore would support the manuscript's publication in Nature Communications after the following comments have been considered.

Comments
]. I was wondering about the interpretation of the black line shown in Fig. ]b, which shows a pronounced temperature dependence as indicated in Suppl. Fig. Z. I missed assignment of this peak in the text, even though this is the range with the most pronounced T-dependence (at least in the log-x representation).
Presumably the confusion has arisen because of the figure caption. We have added a line to the figure caption explaining that the black line is due to diffusive α relaxation. The main text already contains a paragraph (starting with "The spectra in this low frequency range […]") most of which is devoted to diffusive α relaxation and linking it to the measured temperature-dependent shear viscosity.

Ta. As far as I understand the WAXS and SAXS spectra are "stitched together" just below ~U.T /Å, so it would be good to specify the exact q-range in the methods section.
This is a very good point. The data were stitched together using a single scaling factor for all the SAXS data to produce a good match in the range q = M."NO to M."PH Å -" (i.e., just outside of the range of the boson peak feature). This is now described in some detail in the Methods section.
Tb. I am a bit concerned about the influence of the background either from the capillary or the direct beam. It is mentioned that plastic capillaries were used, so it would be useful to plot the SAXS/WAXS background from the capillaries (without any sample) to ensure that there is no matching peak at this momentum transfer range.
The graph below shows the background compared to the signal (and fit) at low temperature showing the bosonpeak feature. The amorphous polycarbonate capillary has a relatively low scattering cross section. Most importantly, it is flat in the region around M.HN Å -" . The sharp peak at M.HU Å -" is due to scattering from the mica windows on the Linkam stage, which remain at room temperature, and is therefore not seen to change with temperature. We are therefore confident that this feature is not caused by scattering from the capillary.
Tc. Suppl. Fig. f needs some additional care. It would be helpful to make the figure labels visible, include proper legends and also re-bin (using larger q-bins) the SAXS data to improve the readability of the data in the q-range above U.U]/Å.
Supplementary Figure P has been improved with bigger and more informative labels, and a larger size. Because of the large amount of data ("MZ temperatures for each WAXS and SAXS), we have not been able to re-bin (it appears that DAWN only allows re-binning of one data set at a time) or re-colour in a better way. Of course, the raw data will be made available in a public data repository on publication. Na. While indeed as it is pointed out in the paper, the pre-peak is reproduced in the Si-Si S(q), similar peaks in the region q<]/Å with even stronger temperature dependence are shown in the S(q)'s of C-O, C-Si and also in the O-O. Do these temperature-dependent signals reflect changes intermolecular or intramolecular distances? I.e. it should be specified whether in the calculation is included atoms within the same molecule or only intermolecular distances. 92º-110ºC averaged, ×10 t, ×10 background Z Indeed, by definition, the Si-Si S(q) exclusively contains information about the inter-molecular distances (as we have one Si atom / TBOS molecule). As such, it is reasonable to expect the emergence of similar pre-peaks in the other S(q) as well, given that we did include both intra-and inter-molecular contributions in the computation of the S(q). To highlight this aspect, we have now computed, for the O-O case, an S(q) obtained from the PCF relative to the O-O inter-molecular contributions only. It can be seen from Figure R" (which we have also added to the Supplementary Information) that the S(q) (relative to the O-O intermolecular contributions only) resembles (as it should) the Si-Si S(q) and similarly displays the emergence of a pre-peak. We have also specified in the revised version of the manuscript that our PCFs and thus S(q) have all been obtained by considering all (i.e., inter-as well as intra-molecular) contributions. Nb. It might be helpful to plot all these results in Suppl. Fig. ]X with logarithmic x-axis (as in Fig. Xb).
We agree and have updated Supplementary Figure "b (now "O) in the revised version of the manuscript as suggested.
Nc. The peak at q~U.T /Å is attributed to clusters with size ~Tπ/q = TZ.V Å. Would then one not expect a similar peak in the Si-Si PCF at this range (Suppl. Fig.]

j)? Or is this range limited by the simulation box size?
Thank you for raising this important point. The size of the simulation box (~dN Å at UM K) is sufficient for us to observe the emergence of these clusters. As illustrated in Figure RH(a), these structural features emerge between M.H and M.b Å -" in our simulated Si-Si S(q) (here, we focus on Si alone as its structural signatures are the most representative of the topology of the TBOS network), which translates to a real-space range between "b and HN Å. The signature in the Si-Si PCF is present-the structure factor is basically the Fourier transform of the PCF after all-and it is highlighted in Figure RH (b). To be sure that the emergence of the pre-peak in our simulations is genuine: • We have thoroughly investigated the influence of sampling (i.e., how many frames we considered within a certain time interval) and binning (i.e., the "resolution" of the PCF), and found our results to be "MM% robust in that respect. • We have invested a significant amount of computational time to generate three additional, statistically independent models of TBOS, following the same computational protocol discussed in the main text. The total as well as the Si-Si S(q), computed as a function of temperature for the four models is reported in Figure RH(a). Note that, in addition to the (expected) variability of the S(q) as a whole (on our scale we are bound to generate/observe an ensemble of slightly different amorphous systems) the pre-peak is present in all the four models. As such, we are now entirely confident in claiming that the emergence of the pre-peak is indeed captured by our molecular dynamics simulations. We have now added these data into the Supplementary Information and we have highlighted the results of the additional models in the main text as well.

X. Concerning the Voronoi analysis:
Xa. It would be helpful to include some additional explanations concerning the formalism used in the corresponding indices and how one can derive the corresponding symmetry.
Thank you for bringing this matter to our attention. In revising our results, we realised that, in the original version of the manuscript, we inaccurately labelled the indices of some of the Voronoi polyhedra. This is because, upon removing the small faces (i.e., faces with area < "% of the total area of the Voronoi polyhedron in question) from our Voronoi polyhedra, voro++ (the library we originally used for this calculation) does not update the number of edges and/or vertices for each face. As the Schlafli notation we have used to identify the Voronoi polyhedra consists in a tuple of eight indices (<n#, n$, n%, n&, n', n(, n), n"*> as discussed in the Supplementary Note Z), where ni refers to the number of faces with i vertexes/edges, the labelling of the Voronoi polyhedra was not guaranteed to be correct.
The solution to this issue is to exclude the neighbouring particles responsible for the emergence of the small faces and to re-compute the Voronoi polyhedra in their absence. We have now implemented this strategy, and we have also cross-checked the robustness of our results by utilising a second piece of software, Ovito. " As such, we have updated Figure N(c) to reflect the new/correct Voronoi indices.
We stress that the only result that has been impacted by this pitfall is the Voronoi indices: every other feature of the Voronoi polyhedra remains unchanged.
In terms of the relationship between the Voronoi indices and the symmetry of the corresponding Voronoi polyhedra, we have now compiled a rather comprehensive list, summarised in the table below (which has also been added to the revised version of the Supplementary Information). In addition to the corresponding coordination number (CN) for each Voronoi index tuple, we have grouped the different Voronoi polyhedra according to the shortrange order they have been reported to be indicative of. Xb. It would be interesting to include a schematic of the Si-SI network and polyhedral derived from the simulations to facilitate the interpretation.
We have tried to follow this suggestion, but the very disordered nature of the Si-Si network within TBOS is such that any meaningful attempt to visualise its topology results in a rather messy/uninformative picture (see Figure  RZ). We hope that the additional information/discussion about the short-range order of the Si-Si network highlighted in the next point (point bc -see the answers to Reviewer #H as well, for additional structural/orientation analysis) would serve to provide a clearer picture of the complex topology of this system.

Xc. The analysis indicates that network of the Si atoms exhibits icosahedral coordination at any given temperature, whereas the observed truncated octahedra emerge near the T_g. One of the main conclusions of the paper is however that the observed transient clusters can inhibit crystallization. According to this interpretation, wouldn't one expect that the cluster formation should emerge in the supercooled regime?
In light of the revised indices for some of our Voronoi polyhedra (see point ba), we had to reconsider our original claim regarding the emergence of icosahedral and BCC order. As can be seen from Figure Rb polyhedra as well, whose population decreases as we approach the Tg. On the other hand, the fraction of (defective) HCP polyhedra (<M,b,b,H>) appears to increase just below Tg, albeit it sharply decreases to a negligible proportion at ~"HM K. We note that the emergence, below Tg, of the over-coordinated molecules (<M,b,b,Z,Z>, <M,b,d,b,H> and <M,b,d,H,H,">) we observed in the original version of the manuscript, persists even in light of our revised analysis. Whilst it is true that, even in metallic glasses (characterised by a much simpler topology if compared to TBOS), the fraction of crystalline-like Voronoi polyhedra is usually minimal ("-N %, see e.g. Ref. 15 ), the variety of Voronoi O polyhedra we observed for TBOS (~"M # ) is such that we can only rely on averaged information, such as the number of faces and the volume of the polyhedra (see Figure Na and Nb), to investigate statistically significant trends in terms of the topology of TBOS as a function of temperature.
To this end, we have leveraged the framework of Ref. "# to investigate the local molecular environments of the Si network within our TBOS models. The results are summarised in Figure Rb (b), and clearly exclude the emergence of either FCC or BCC structures at any temperature. Conversely, a small fraction (~Z%) of TBOS molecules sit in HCP environments, whose abundance decreases as the system is cooled below the glass transition.
Further evidence for the absence of substantial crystalline order is provided by the bond order analysis illustrated in Figure Rb(c), where we monitor the values of the Steinhardt parameters q$ and q& as a function of temperatures (again, for each Si atom within the TBOS network). We remark that we have employed the local average version of q$ and q& pioneered in Ref. "$ . It can be seen that the system is characterised by values of q$ that are characteristic of HCP order. However, the fact that the values of q& are consistently within the liquid regime is further proof of the fact that this tendency toward the HCP structure is entirely frustrated by the emergence of over-coordinated environments, which add even more disorder into the network. Note that the supercooled liquid is in fact more consistent with HCP order (in terms of q$) than the glassy state, which is entirely consistent with the local environment analysis of Figure Rb

Reviewer #2 (Remarks to the Author):
The work by González-Jiménez studied the dynamical properties, i.e. boson peak, of a molecular glass by combining experiments and simulations. A link to the local structure is suggested. However, this paper is more technical relevant with lots of assumptions involved in the analysis, so it is not convincing to me. Therefore, I cannot recommend its publication in Nat. Commun., a more specialized journal may be more suitable. See some detailed comments below.
]) The introduction should be carefully re-organized to target the problem to be solved in this paper. For example, the main results are related to boson peak and glass structure, but the introduction mainly focused on different relaxation modes, like beta relaxation. Not many details about the boson peak are shown.
We argue in the introduction that both secondary relaxations (Johari-Goldstein ones at least) and the boson peak are related to intermolecular interactions and therefore the intermolecular structure and dynamics. In our view that merits the one paragraph on secondary relaxations in the introduction. (The absence of) β relaxations in the OKE spectra are discussed again in the Discussion and conclusions, so it is important to keep that introductory paragraph.
To strengthen the paragraph on the boson peak, we have added references to the recent paper on the origin of the boson peak in Nature Physics (HMHH), which is also discussed briefly in the conclusion section. In addition, we have added a reference to the recent (HMHH) book "Low-Temperature Thermal and Vibrational Properties of Disordered Solids" by Miguel A. Ramos. Finally, a couple of sentences have been added regarding the VDOS in glasses vs. crystals.

T) In the introduction, "Below the glass-transition, secondary or β-relaxation processes occur due to faster dynamic processes" may not be necessarily true. I believe different relaxation modes should coexist always in a thermodynamic state, but their observation will depend on the experimental time scale.
We agree on this but feel that a discussion on observation timescales is not appropriate for this manuscript on the boson peak. However, the word "occur" does indeed suggest that β relaxations only come into existence below Tg, which is not necessarily correct. Therefore, we have modified the phrasing in that paragraph.

N) There is a contradiction about the origin of the secondary relaxations in the introduction, whether intramolecular or intermolecular.
P There is no contradiction: β relaxations can be intermolecular or intramolecular and which one is dominant depends on the system. In polymers, intramolecular (small angle orientational diffusion of the side chains or functional groups) β relaxations are thought to be dominant whereas in small-molecule and metallic glasses intermolecular/atomic β relaxations (or Johari-Goldstein) are dominant. Since this manuscript is on a smallmolecule glass former, the discussion in this paragraph is biased towards Johari-Goldstein β relaxation.
X) The detailed data analysis from the experimental spectra depends on the assumption that they will be simplified for some molecules with ordinary symmetry, tetrahedral and octahedral. This has not been rationalized. Also, in the glass state, there are enormous types of these local structures with different distortions, how would this fact influence the results?
The fact that the depolarised Raman spectrum depends on the anisotropic part of the molecular polarisability tensor is in basic textbooks on Raman scattering while the vanishing of the anisotropic part in molecules with tetrahedral symmetry is basic group theory. However, it is of course correct that a molecule such as TBOS will be distorted and will not have perfect tetrahedral symmetry. The effect of this distortion is discussed in detail on page Z (plus more details in the Methods section) where we discuss DFT calculations of the molecular polarisability tensor for Z",NMM low-energy conformers.
j) The spectrum fitting with four functions of many free parameters seems highly possible to bring over-fitting. This is very critical for the overall interpretation of the results. This is of course correct in principle. However, TBOS was specifically chosen as its near-tetrahedral symmetry simplifies the spectra. As a result, the spectra taken in the glassy state ( Figure "(c)) clearly show Z peaks below b THz, making fitting very straightforward and unambiguous. Well above the glass transition, these very same Z peaks can still be readily observed, see, for example, Supplementary Figure N, ZMM K. At still higher temperatures, the diffusive α relaxation band shifts to ever higher frequencies, making fitting trickier. However, as Supplementary Figure O shows, the fitted relaxation time-constant closely follows the Stokes-Einstein law using the viscosity measured with rheometry as expected. Therefore, one can have full confidence in these fits.

Z) The structure analysis by pair-correlation function and the Voronoi tessellation does not account for the orientation information in the local structures, which should be important in this molecular glass.
TBOS molecules can be considered as approximately tetrahedral units. The first question in terms of the orientation of these tetrahedral units is whether, on average, they would assume any preferential orientation with respect to a given reference axis, particularly in the glassy state. To investigate this possibility, we have now computed, for each TBOS molecule, the angle Θz formed by each of the four Si-O molecular axes (defined as the vectors originating from the Si atom and pointing toward each of the four oxygen atoms) and the Cartesian z-axis. We have then averaged these four angles over all the molecules within a given frame to obtain an average orientation <Θ>z. The probability density function of <Θ>z, computed over "MMM statistically independent frames within "M ns long trajectories, is reported in Figure RN(a) as a function of temperature. It can be seen that the distributions peak at ~UM o at every temperature, which demonstrates that the system is isotropic even in the glassy state. Note the narrowing of the distribution as the temperature decreases -which is indicative of the much slower rotational dynamics within the glassy state.
Next, we wanted to investigate the relative orientation of the TBOS tetrahedral units within their local environments. To this end, we define, for each molecule, the plane illustrated in Figure Rb (b). We then leverage the so-called Solid Molecule Angle Criteria (SMAC) order parameter, described in Ref. "& . In a nutshell, SMAC measure the degree to which molecules, within a certain distance from each other, are similarly oriented. In this case, we considered up to the H nd coordination shell of the TBOS molecules, which encompasses ~"H neighbours. The orientation of each molecule is defined by the vector normal to the plane illustrated in Figure Rb (b). Higher (lower) values of SMAC correspond to more (less) ordered molecular environment in terms of the relative orientation of, in this case, the tetrahedral SMAC units.
The results are summarised in Figure Rb(c), where we have averaged the value of SMAC for each molecule within a given frame (which yields <SMAC>) and then built a probability density function over different frames (same sampling as Figure Rb(a)). Interestingly, it appears that the glassy state is less ordered, in terms of the relative orientation of the TBOS molecules, than the supercooled liquid. This result is consistent with the progressive loss U of short-range order we have observed via the analysis reported in Figure RZ, and provides further evidence of the fact that the over-coordinated (Voronoi polyhedra characterised by > "H faces, particularly "N and "d faces) clusters of TBOS molecules we have observed in the glassy state are not characterised by any degree of crystalline order. This analysis has now been added to the Supplementary Information. p) The direct comparison of the simulations by a simplified model with the experiments is not convincing. The consistency of the partial structure factor only shows the overall averaged information similarly, not necessarily mean the model it's correct.
The MD simulations accurately reproduce: (") the temperature-dependent behaviour of the OKE spectra; (H) the temperature-dependent behaviour of the WAXS data, including the emergence of the pre-peak at ~M.H Å -" ; (Z) the temperature-dependent changes of the relative population of gauche/trans configurations in the alkoxide chains of TBOS-an additional piece of analysis we have now included in the Supplementary Information, see our answer to Reviewer #Z. As such, the agreement between simulations and experiments is remarkable, especially if we consider that the force field we have utilised was not explicitly parametrised with respect to alkoxides such as TBOS. In addition to this, we have now invested a substantial amount of computational time to generate three additional, statistically independent models of TBOS, which allowed us to strengthen even further our claims in terms of the structural features we observe in our molecular dynamics simulations (see the answer to Reviewer #", point Zc).

f) At elevated temperatures, anharmonicity is important, which should have been into consideration in the spectrum analysis.
Our simulated OKE spectra, which-because of the inherent symmetry of the TBOS molecules are directly related to the translational motion of Si atoms-have been obtained from the Fourier transform of the Si-Si velocity autocorrelation functions, do take into account anharmonicity, given that our MD simulations have been performed at finite temperature.

Reviewer #3 (Remarks to the Author):
The manuscript of González-Jiménez et al. entitled "Understanding the emergence of the boson peak in molecular glasses" aims at unveiling the boson peak in a highly symmetric glass former using femtosecond optical Kerr-effect spectroscopy as the central characterization technique, at different temperatures. They then use other experimental techniques (WAXS, SAXS, Raman) and simulation to characterize this boson peak. It is because the characterized molecule, TBOS, is symmetric that the authors can tackle this longstanding problem of the origin of the boson peak.
Two intermolecular bands in the OKE spectrum have been mainly analyzed. In Figure  ] d, the amplitudes of both bands are shown versus the temperature, and in Figure T a, their ratio is displayed. Through theses analyses, the authors assign the low frequency band as the boson peak, and the high frequency band as the acoustic phonon band. Thanks to this assignment, they surmised that the boson peak is associated with the formation of clusters of TBOS molecules undergoing collective oscillation. WAXS experiments lead to a distance of TZ.V Angstrom for clusters, and from MD simulation, Tj.] Angstroms. Conclusions drawn from the OKE experiments can have been made because of the symmetric nature of the TBOS molecule. Regarding the over-coordinated molecules in the clusters, deduced from MD, the authors wrote: "the MD simulations are strongly suggestive that these structures are clusters of over-coordinated TBOS molecules". For an article in Nature Communications, the authors must be more convincing. For instance, as the molecular behavior can be probed, the over-coordination can be proved.
We have in fact conclusive evidence, from our molecular dynamics simulations, of the emergence of overcoordinated TBOS molecules -as proven by the Voronoi analysis reported in Fig. Na. We also have conclusive evidence for the fact that these over-coordinated molecules tend to cluster in aggregates whose size increases as the system is cooled into the glassy state. This was mentioned, briefly in the original version of the manuscript, but we have now included the probability density function of the largest cluster (consisting of TBOS molecules characterised by Voronoi polyhedra with "N or more faces) in each frame of our MD trajectories in the Supplementary Information. The result is summarised in Figure Rd as well, as a function of temperature. Finally, we have shown ( Fig. N (d)) that over-coordinated TBOS molecules (particularly those characterised by Voronoi polyhedra with "d faces) are responsible for populating the high-frequency region of the Fourier Transform of the velocity-velocity autocorrelation function of the Si atoms (a quantity which is directly related to the experimental OKE spectrum given the symmetry of the TBOS molecules). In light of all of the above, we are quite confident that these particular structural features do correspond to those identified by the OKE measurements. As such, we have replaced the rather cautious wording, "strongly suggestive" with "show" in the revised version of the manuscript. Our apologies! In the original manuscript the computed spectrum was simply shifted and therefore appears (if you look at the peak wavenumbers rather than the shape) to be completely mismatched. It is common for normal mode calculations to slightly overestimate the normal-mode frequencies. Therefore, we have matched the calculated frequencies of the CH-stretch modes with the experimental ones by multiplying them with M.UZd. This shifts the calculated Si-O-C modes down a bit as well and the correspondence is now much more obvious. The figure and its caption have been updated.
The figures below show the spectra from dON cm -" (the lowest we have been able to go experimentally with the available IR windows) to Z,MMM cm -" . It can be seen that the correspondence is reasonable over the entire range. The manuscript only shows the ",MMM-",HMM cm -" region as this should show a prominent peak in highly symmetrically coordinated (tetrahedral, octahedral) alkoxides, as it does in TBOS.

Indicate what are qN and qj.
We apologise that several parameters were not defined in the manuscript! Changes were made to Supplementary note H (Analysis of WAXS data) and the captions/headers of Supplementary Tables H and Z. All the amplitude (Li and Gi), width (γi and σi), and position (qi) parameters are now properly defined and the definitions referred to in the text.

Below ]]U K: the value seems to be ]]U K.
Agreed, changed to "at".

It is written that "the TBOS molecules reduce the number of gauche defects in the alkoxide side chains". Such a behaviour should be observed in MD.
Indeed, this behaviour is captured by the MD simulations. To illustrate this aspect, we have computed the probability density function of the value of all the H-C-C-H dihedral angles ( ) within the alkoxide side chains, averaged across bM ns at each temperature. The results are summarised in Figure RP. As the temperature decreases, the distribution of narrows (as the mobility of the side chains is greatly reduced in the glassy state) and we observe a reduction of the relative population of gauche configurations with respect to trans configurations (as demonstrated by the ratio between the height of the PDF peaks at ~"HM° [trans] and ~dM° [gauche], which increases from ~H.N at bbM K to ~Z.b at UM K). Trans configurations are clearly more energetically favoured at any temperature, potentially due to steric effects. Figure RP has now been added to the Supplementary Information. p.j: The authors consider the two coordination shells in the RDF, computing T and ]T molecules respectively. According to these spectra, these numbers should change with the temperature. Some relevant results can be obtained.
Indeed, the Si-Si coordination number of the system does change with temperature, as illustrated in Fig. bd. One aspect we failed to stress in the original version of the manuscript is that the short-/medium-range order in supercooled liquid, and even more so in glassy TBOS, is rather peculiar. In fact, we observe a first Si-Si coordination shell (up to O.N Å) that on average includes only ~H.N nearest neighbours. This is interesting, as even in the glassy state the TBOS molecules can almost never quite arrange themselves into a b-coordinated configuration within their first coordination shell. Instead, they either assume a H-fold or at most a Z-fold coordinated configuration, illustrated in Figure RU, where the TBOS molecules are almost lying in the same plane. We argue that these peculiar structural features have to do with steric effects related to the presence of the alkoxide chains (on top of the fact that tetrahedra do not fill ZD space), which prevents a denser packing of the system. Rather distinct from this under-coordinated short-range feature, we also observe a second Si-Si coordination shell (which extends up to ~"H Å) that includes, on average, ~"H nearest neighbours (consistently with our Voronoi analysis, see e.g. Fig. Na). We have now highlighted this structural feature of the system in the revised version of the manuscript. When the authors consider random positions within a cubic box: what is the random algorithm used? It is necessary to specify that the periodic boundary conditions have been applied.
We have used Packmol 17 to generate the initial random arrangement of TBOS molecules within the simulation box. Details about the algorithm used by this particular piece of software can be found in Ref. 18 . We have added this information in the Methods section. However, we note that, given we have thoroughly equilibrated the system in a highly diffusive regime, the initial (random) configuration of the TBOS molecules within the simulation box does not have any impact on our subsequent molecular dynamics simulations.
Simulations are performed in the NPT ensemble, meaning that a thermostat and a barostat are needed simultaneously. So, why the authors wrote: "to sample the canonical ensemble"?
Thank you for pointing this out. We have now changed this sentence: "The Bussi-Donadio-Parrinello thermostat #$ was used to sample the canonical ensemble, in conjunction with the Berendsen barostat. %& " to: "The Bussi-Donadio-Parrinello thermostat #$ was used in conjunction with the Berendsen barostat %& to sample the isobaric-isothermal (NPT) ensemble." The authors wrote: "In this regime, the self-diffusion is sufficient to sample a reasonable portion of the configurational space for the liquid phase". One point in the configurational space corresponds to fp,jjT coordinates. What do the authors mean with this reasonable portion?
We apologise for the confusion. We have computed the self-diffusion coefficient of the Si atoms at that temperature to be ~" nm 0 /ns at bbM K, which indicates that the system is in a highly diffusive, hydrodynamic regime (we have verified this by analysing both the mean square displacement and the self part of the intermediate scattering function, which at that temperature displays a single exponential decay characteristic of a purely hydrodynamic regime). As we have equilibrated the system for > HM ns at that temperature, we are confident that we have generated statistically independent configuration of the liquid to serve as starting points for our quench/cooling protocol.
To make this clearer, we have now changed the sentence: "In this regime, the self-diffusion is sufficient to sample a reasonable portion of the configurational space for the liquid phase, thus providing an adequate starting point for the cooling ramp." to: "Under these conditions, the self-diffusion coefficient of the Si atoms was computed to be ~" nm 0 /ns indicative of a highly diffusive hydrodynamic regime. As the system was equilibrated for a further HM ns at that temperature and pressure, statistically independent well-equilibrated configurations of the liquid were generated providing an adequate starting point for the cooling ramp." in the revised version of the manuscript.
Finally, how can the study presented in this paper be correlated with the recent paper in Nature Physics on the origin of the boson peak?
Thank you. Of course, this nice paper only came out when this manuscript was submitted. A sentence has been added to the Discussion linking our results to the recent theoretical work by Hu and Tanaka in Nat. Phys.